function [V, T, TE, E, ET] = refine_mesh_uniform(V, T, TE, E, ET)
% function [V, T, TE, E, ET] = refine_mesh_uniform(V, T, TE, E, ET)
%
% this version can update the information on E automatically!
% which save large amount of time costs
%
n_elem = size(T  ,1);
n_vert  = size(V  ,1);
n_edge = size(ET,1);
% 1: the father edge;  2: the local neighboured triangle;  

% all new vertices:  
V = [V; (V(E(:,1),:) + V(E(:,2),:))/2];
    
% pick up vertices' index element-wisely
v1 = T(:,1); v2 = T(:,2); v3 = T(:,3);
v4 = n_vert + TE(:,1);
v5 = n_vert + TE(:,2);
v6 = n_vert + TE(:,3);

% all new edges, 2*n_edge + 3*n_elem
vc = n_vert + (1:n_edge)';
ET = E;  % just save the old edge for temperal using
         % note that ET is generated only according to TE at last!
E = [E(:,1),  vc; ...   % first two column in edge order
        E(:,2),  vc; ...
        v5     ,  v6; ...  % the rest inner edge is in element order
        v4     ,  v6; ...
        v4     ,  v5];

% pick up new subedge's position element-wisely
% :::Remember that the index of the old Edge (now stored in ET) is chaos!
E1 = TE(:,1); E2 = E1 + n_edge;   % by default is anti-clock wise
idx = (ET(E1,1) == T(:,3));       % find the abnormal
tmp = E1(idx); E1(idx) = E2(idx); E2(idx) = tmp;  % and exchange
E3 = TE(:,2); E4 = E3 + n_edge;    % by default is anti-clock wise
idx = (ET(E3,1) == T(:,1));       % find the abnormal
tmp = E3(idx); E3(idx) = E4(idx); E4(idx) = tmp;  % and exchange
E5 = TE(:,3); E6 = E5 + n_edge;    % by default is anti-clock wise
idx = (ET(E5,1) == T(:,2));       % find the abnormal
tmp = E5(idx); E5(idx) = E6(idx); E6(idx) = tmp;  % and exchange
% 

% all new triangles
T1 = [v1, v6, v5];
T2 = [v2, v4, v6];
T3 = [v3, v5, v4];
T4 = [v4, v5, v6];
T   = [T1; T2; T3; T4];


E7 = 2*n_edge + (1:n_elem)';
E8 = E7 + n_elem;
E9 = E8 + n_elem;

% new TE informations
TE1 = [E7, E4, E5];
TE2 = [E8, E6, E1];
TE3 = [E9, E2, E3];
TE4 = [E7, E8, E9];
TE = [TE1; TE2; TE3; TE4];

% new ET informations, 
n_new_edge = n_edge*2 + n_elem*3;
n_new_elem = n_elem*4;
ET = zeros(n_new_edge, 2);
idx =  ones(n_new_edge, 1);
for tri = 1: n_new_elem
    bnd1 = TE(tri,1);
    ET(bnd1, idx(bnd1)) = tri;
    idx(bnd1) = idx(bnd1) + 1;  % do not set to 2 directly
    bnd2 = TE(tri,2);
    ET(bnd2, idx(bnd2)) = tri;
    idx(bnd2) = idx(bnd2) + 1;  % do not set to 2 directly
    bnd3 = TE(tri,3);
    ET(bnd3, idx(bnd3)) = tri;
    idx(bnd3) = idx(bnd3) + 1;  % do not set to 2 directly
end


% function [node, elem] = refine_mesh_uniform(node, elem, n)
% % function [node, elem] = refine_mesh_uniform(node, elem, n)
% % This function come from AFEM by LongChen
% % UNIFORMREFINE refines the current triangulation by dividing
% %   each triangle into four similar triangles
% % 
% 
% %fprintf('\n Maping solution: current level: level %g  to level %g ...',old_level,new_level);
% %--------------------------------------------------------------------------
% % Construct data structure
% %--------------------------------------------------------------------------
% for i=1:n
% 
% edge = [elem(:,[1,2]); elem(:,[1,3]); elem(:,[2,3])];
% edge = unique(sort(edge,2),'rows');
% N = size(node,1); NT = size(elem,1); NE = size(edge,1);
% d2p = sparse(edge(:,[1,2]),edge(:,[2,1]),[1:NE,1:NE],N,N);
% 
% %--------------------------------------------------------------------------
% % New nodes from the mid points of each edge
% %--------------------------------------------------------------------------
% 
% 
%     
% 
% newnode = (node(edge(:,1),:)+node(edge(:,2),:))/2; 
% node = [node; newnode]; 
% 
% marker = N+1:N+NE; 
% 
% %--------------------------------------------------------------------------
% % refine each triangle into four triangles
% %     3
% %    / \
% %   6 - 5
% %  / \ / \
% % 1 - 4 - 2
% %--------------------------------------------------------------------------
% for t=1:NT
%     p(1:3) = elem(t,[1 2 3]); 
%     p(4) = marker(d2p(elem(t,1),elem(t,2)));
%     p(5) = marker(d2p(elem(t,2),elem(t,3)));
%     p(6) = marker(d2p(elem(t,3),elem(t,1)));
%     elem(t,[1 2 3]) = [p(1) p(4) p(6)];
%    elem(size(elem,1)+1,[1 2 3]) = [p(4) p(2) p(5)];
%     elem(size(elem,1)+1,[1 2 3]) = [p(6) p(5) p(3)];
%    elem(size(elem,1)+1,[1 2 3]) = [p(4) p(5) p(6)];
% end
% 
% 
% end
%-------------------------------------------